﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;




using HA.ModFlow.Domain.DomainLayer;
using HA.ModFlow.Infrastructure;





namespace HA.MT3DMS.DomainLayer
{




    public static   class ComputeMass
    {

        #region 字段属性


        #endregion



        #region 共有方法

        /// <summary>
        /// 非用户选定计算单元的模式  根据  IBound  计算所有网格
        /// </summary>
        /// <param name="xLength">xLength和yLength分别按坐标从小到大存放网格坐标系下的网格边框点的坐标（存放坐标不合适，应为可能旋转平移了）</param>
        /// <param name="yLength"></param>
        /// <param name="zLength">zLength中的值从上到下存放网格中心点的坐标 </param>
        /// <param name="porosity"></param>
        /// <param name="calMin"></param>
        /// <param name="calMax"></param>
        /// <param name="dsHead">给定的应力期的modflow输出的水头文件数据</param>
        /// <param name="dsConc">给定的应力期的mt3d输出的污染物的timeStep的结果数据</param>
        /// <param name="factory"></param>
        /// <returns></returns>
        public static double Compute(
            double[] xLength,
            double[] yLength,
            List<double[,]> zLength,



            List<double[,]> porosity,

            double calMin,
            double calMax,

            CellCenterDataSet3d dsHead,
            CellCenterDataSet3d dsConc,

             double factory

            )
        {

            #region 进行一些必要的数据有效性的检验


            if (xLength == null || yLength == null || zLength == null || dsHead == null || porosity == null || dsConc == null)
            {
                return 0;
            }

            if (calMin > calMax)
            {
                return 0;
            }

            if (porosity.Count <= 0 || zLength.Count <= 0)
            {
                return 0;
            }


            #endregion


            try
            {

                string timeStepLabel = dsConc.Label;// 标签   格式： MT3D: 应力期_时间步长_传输步长_时间   ModFlow:应力期_时间步长

                double timeStepValue = 0.0;



                #region 求得 timeStepValue  该部分能够和下面的方法中的对应部分合并 减少代码的冗余  呵呵

                List<int[,]> iBound = dsConc.EffIBound;//取得mt3d污染物结果中的单元格的有效性标志

                for (int layer = 0; layer < iBound.Count; layer++)//遍历所有单元格
                {
                    for (int row = 0; row < iBound[layer].GetLength(0); row++)
                    {
                        for (int col = 0; col < iBound[layer].GetLength(1); col++)
                        {

                            #region

                            double calculationValue = dsConc.Datas[row, col, layer];

                            if (iBound[layer][row, col] != 0 &&  Math .Abs (calculationValue - dsConc.InvalidValue)>1e-8 )//单元格有效并且 是 有效值
                            {

                                if (calMin <= calculationValue && calculationValue <= calMax)//污染物结果文件中的数据是否在指定的  CalculationValueRange 范围内
                                {


                                    //单元格体积
                                    double cellVolume = (yLength[row + 1] - yLength[row]) * (xLength[col + 1] - xLength[col]) * (zLength[layer][row, col] - zLength[layer + 1][row, col]);//单元格体积



                                    double iCellHeadValue = dsHead.Datas[row, col, layer];


                                    //记录该单元格的结果 Mass
                                    double cellMass = 0;

                                    if (iCellHeadValue < zLength[layer][row, col])//水头值<顶部高程，需要乘以比例
                                    {
                                        cellMass = cellVolume * porosity[layer][row, col] * calculationValue * factory;
                                        cellMass = cellMass * (zLength[layer][row, col] - iCellHeadValue) / (double)(zLength[layer][row, col] - zLength[layer + 1][row, col]);

                                    }
                                    else
                                    {

                                        cellMass = cellVolume * porosity[layer][row, col] * calculationValue * factory;

                                    }

                                    //累加该应力期下所有单元格的值

                                    timeStepValue = timeStepValue + cellMass;

                                }


                            }// if (iBound[layer][row, col] != 0)//单元有效

                            #endregion



                        }
                    }



                }


                #endregion


                return timeStepValue;

            }
            catch
            {

                return 0;
            }


      

        
        }


       /// <summary>
       /// 用户指定计算单元格
       /// </summary>
        /// <param name="xLength">xLength和yLength分别按坐标从小到大存放网格坐标系下的网格边框点的坐标（存放坐标不合适，应为可能旋转平移了）</param>
        /// <param name="yLength"></param>
        /// <param name="zLength">zLength中的值从上到下存放网格中心点的坐标 </param>
        /// <param name="porosity"></param>
        /// <param name="calMin"></param>
        /// <param name="calMax"></param>
        /// <param name="dsHead">给定的应力期的modflow输出的水头文件数据</param>
        /// <param name="dsConc">给定的应力期的mt3d输出的污染物的timeStep的结果数据</param>
        /// <param name="factory"></param>
       /// <param name="userSpecifiedCellIJK">用户指定的计算的单元格  数组中存放顺序为 Row，Col，Layer</param>
       /// <returns></returns>
        public static double Compute(
            double[] xLength,
            double[] yLength,
            List<double[,]> zLength,


            List<double[,]> porosity,

            double calMin,
            double calMax,

            CellCenterDataSet3d dsHead,
            CellCenterDataSet3d dsConc,

             double factory,

            List<int[]> userSpecifiedCellIJK

            )
        {

            #region 进行一些必要的数据有效性的检验



            if (xLength == null || yLength == null || zLength == null || dsHead == null || porosity == null || dsConc == null || userSpecifiedCellIJK == null)
            {
                return 0;
            }
            if (calMin > calMax)
            {
                return 0;
            }


            if (porosity.Count <= 0 || zLength.Count <= 0 || userSpecifiedCellIJK.Count <= 0)
            {
                return 0;
            }


            #endregion

            try
            {

                string timeStepLabel = dsConc.Label;// 标签   格式： MT3D: 应力期_时间步长_传输步长_时间   ModFlow:应力期_时间步长

                double timeStepValue = 0.0;



                #region


                List<int[,]> iBound = dsConc.EffIBound;


                for (int j = 0; j < userSpecifiedCellIJK.Count; j++)//遍历用户选择计算的单元格
                {
                    int row = userSpecifiedCellIJK[j][0];
                    int col = userSpecifiedCellIJK[j][1];
                    int layer = userSpecifiedCellIJK[j][2];




                    #region

                    double calculationValue = dsConc.Datas[row, col, layer];

                    if (iBound[layer][row, col] != 0 && Math.Abs(calculationValue - dsConc.InvalidValue) > 1e-8)//单元格有效并且 是 有效值
                    {

                        if (calMin <= calculationValue && calculationValue <= calMax)//污染物结果文件中的数据是否在指定的  CalculationValueRange 范围内
                        {


                            //单元格体积
                            double cellVolume = (yLength[row + 1] - yLength[row]) * (xLength[col + 1] - xLength[col]) * (zLength[layer][row, col] - zLength[layer + 1][row, col]);//单元格体积



                            double iCellHeadValue = dsHead.Datas[row, col, layer];



                            //记录该单元格的结果 Mass
                            double cellMass = 0;

                            if (iCellHeadValue < zLength[layer][row, col])//水头值<顶部高程，需要乘以比例
                            {
                                cellMass = cellVolume * porosity[layer][row, col] * calculationValue * factory;
                                cellMass = cellMass * (zLength[layer][row, col] - iCellHeadValue) / (double)(zLength[layer][row, col] - zLength[layer + 1][row, col]);

                            }
                            else
                            {

                                cellMass = cellVolume * porosity[layer][row, col] * calculationValue * factory;


                            }

                            //累加该应力期下所有单元格的值

                            timeStepValue = timeStepValue + cellMass;

                        }


                    }// if (iBound[layer][row, col] != 0)//单元有效

                    #endregion





                }



                #endregion



                return timeStepValue;


            }
            catch
            {

                return 0;
            }


        }




        #endregion



        #region 内部方法




        #endregion





    }
}
